function Mvect = steadySolution(Bx, By, Bz, T1, T2, R, M_inf, gamma)
%STEADYSOLUTION Summary of this function goes here
%   Detailed explanation goes here
    B=norm([Bx, By, Bz])+eps; Omega = gamma*B;
    denormenator = (T2+T1)*(T2+T1+2*R*T1*T2)+T1*T1*T2*T2*(R*R+Omega*Omega);
    MM = M_inf*R*T1/(1+R*T1);

    Mx = Bx*Bz*T1*(T2+T1+T1*T2*(R+Omega^2*T2))...
        /( B*B * denormenator ) * MM;

    My = -gamma*Bx*R*T1*T1*T2*T2 / denormenator * M_inf;

    Mz = ( Bz*Bz* (T1*T1+T1*T2*(1+R*T1))+T1*T2* B*B *(1+R*T1)+ B*B*T2*T2*(1+2*R*T1+T1*T1*(R*R+Bz*Bz*gamma*gamma)) )...
        /( B*B * denormenator ) * MM;
    
    Mvect = [Mx, My, Mz];
end

